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i— i I Abstract 

<n : 

In t/iis work a finite element simulation of the motion of a rigid body in a fluid, with 
■ free surface, is described. A completely general referential description ( of which both 

Lagrangian and Eulerian descriptions are special cases) of an incompressible, New- 
tonian fluid is used. Such a description enables a distorting finite element mesh to 
Q\ ' be used for the deforming fluid domain. 

o\ : 

A new scheme for the linearised approximation of the convective term is proposed 

• | and the improved, second order accuracy of this scheme is proved. A second theo- 

rem, which provides a guideline for the artificial adjustment of the Reynolds number 
^ when applying a continuation technique, is also proved. The most effective means 

of eliminating pressure as a variable and enforcing incompressibility are reviewed. 
' A somewhat novel method to generate finite element meshes automatically about 

^ ■ included rigid bodies, and which involves finite element mappings, is described. 

The approach taken when approximating the free surface, is that it may be treated as 
a material entity, that is, the material derivative of the free surface is assumed zero. 
Euler's equations and conservation of linear momentum are used to determine the 
motion of the rigid body. A predictor-corrector method is used to solve the combined 
sub-problems. The resulting model is tested in the context of a driven cavity flow, 
a driven cavity flow with various, included rigid bodies, a die-swell problem, and a 
Stokes second order wave. 



Keywords: rigid body in a fluid; free surface; incompressible, Newtonian fluid; completely 
general referential description; arbitrary Lagrangian Eulerian; A.L.E.; finite elements. 



2 



S.J. Childs, B.D. Reddy 



1 Introduction 

An enormous diversity of problems subscribe to the basic rigid body-fluid-free surface 
theme (see Figure [I] for a schematic form of the problem of interest). The implications 
of free surfaces and fluid-rigid body interactions are a deforming fluid domain and the 
implementation of the majority of numerical time integration schemes then becomes prob- 
lematic where a purely Eulerian description of the fluid motion has been used. The reason 
is that most numerical time integration schemes require successive function evaluation 
at fixed spatial locations. Meshes rapidly snarl when purely Lagrangian descriptions are 
used. 




Free Surface 



Rigid Body 



Viscous, 

Incompressible 

Fluid 



Figure 1: The Motion of a Rigid Body in a Fluid with Free Surface 

One alternative is to use the completely general referential description proposed in 
Childs 0. Eulerian and Lagrangian references are just two, specific examples of an 
infinite number of configurations over which to define the fields used to describe the dy- 
namics of deforming continua. They are both special cases of a more general referential 
description, a description in which a reference configuration is chosen at will (closely 
related to the so-called arbitrary Lagrangian Eulerian method or A.L.E. method). The 
completely general referential description is implemented in this work. 

Solving for the orientation and position of a rigid body entails solving equations of motion. 
These are derived from the principles of conservation of angular and linear momentum. 
Expressing the principle of conservation of angular momentum in more usable terms is 
not as straight forward as the linear momentum case, and a more appropriate reference 
must be resorted to. This reference is fixed within the rigid body in such a way that 
the origin and centre of mass coincide, as do the axes and principal moments of inertia. 
A complete transformation of all quantities, to this more appropriate reference, is used. 
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This is in contrast to the approach used in deriving the fluid equations, where quantities 
are not perceived in terms of the new reference, but instead merely defined to be functions 
of it. The change of reference is accomplished using a transition matrix and the Coriolis 
theorem. 

The approach taken when approximating the free surface, is that it may be treated as a 
material entity i.e. the material derivative of the free surface is assumed zero. Surface 
tension is not considered significant in the free surfaces investigated. 

A predictor-corrector method is used to solve the combined rigid body, fluid and free 
surface sub-problems. A backward difference scheme is used to approximate the time 
derivative in the fluid sub-problem (in this regard see Childs ||). The finite element 
method is used for the spatial (referential "space") discretisation. Nonlinearity is cir- 
cumvented by way of a new second order accurate linearisation, and the use of Picard 
iteration if necessary. Picard iteration amounts to "linearising" the non-linear term with 
the solution obtained during the previous iterate. A penalty method is used to eliminate 
pressure as a variable and a Q2-P1 element pair is used as a basis. A somewhat novel 
method to generate finite element meshes automatically about included rigid bodies, and 
which involves finite element mappings, is described. 

The more conventional methods for solving initial value problems are used to solve for the 
position and orientation of the rigid body. The definitions of linear and angular velocity 
are coupled to the equations of motion to produce a system of twelve equations (which 
are of course also coupled to the fluid equations) and a Runge-Kutta-Fehlberg method 
is then used. A backward difference scheme and the finite element method are used for 
the respective temporal and spatial discretisations in the free surface sub-problem. 

The overall model was tested in the context of a driven cavity flow, a driven cavity flow 
with various, included rigid bodies, a die-swell problem and a Stokes, second-order wave. 
Results for these examples are presented in Section || The programme was written by 
the author in Fortran and originally run on a model 400 Dec Alpha 3100. 



2 The Equations Governing the Motion of the Rigid 
Body 

What follows is a brief exposition of the equations of motion. This is, of course, review of 
classical work contained in a great variety of references. Its inclusion is deemed appropri- 
ate not for the sake of completeness alone, but because Euler's equations (the equations 
governing the angular motion) and the equations governing the translation are written 
in terms of two distinctly different references. 

The equations of motion, six in all, are essentially an explicit re-write of the principles 
of linear and angular momentum. The equations which govern the translation, 

dVrb _ 

dt m 
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need no introduction, they being a straight forward statement of the principle of conser- 
vation of linear momentum. In the above equation v r b denotes the velocity of the rigid 
body, m denotes the mass, £/ is the sum of the external forces acting on the rigid body 
and t denotes time. Expressing the angular momentum equations in a more practica- 
ble form is not as straight forward and necessitates the change to the reference already 
mentioned. The conservation principle of angular momentum may be stated as 

where q is the angular momentum and S(f A /) is the sum of the torques which arise as 
a result of the external forces acting about the centre of mass (r is the position vector 
of a material point within the rigid body and / is the force). The equations for the 
conservation of angular momentum, though easily understood, are not of much use when 
stated in the above form. 



A More Appropriate Reference 

The fact that 



/ pr dVL r b = 



for a reference whose origin coincides with the centre of mass (Q r b is the rigid body 
domain depicted in Fig. ^) immediately suggests the existence of a more appropriate 
reference than the one in terms of which the above conservation principle is written. For 
such a reference fixed within the rigid body, 



^ = 

dt 



in addition. This suggests expressing 

q = / pr Av dQ r b 

f di *• f * dx 

= / pr A — dQ rb + pr dQ rb A — 

in terms of a more appropriate reference. Converting to this new, more appropriate 
reference by means of the Coriolis theorem, the relation 

q = [ pr A (<jJ A r)dVt rb [ — = for the new reference] 

Jn rb \dt J 

is obtained, where uj is the angular velocity. It is then a fairly elementary exercise to 
show that 



q = Ju), 

TThe "s are resorted to temporarily in anticipation of the change to the new reference. 



(2) 
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where J is the inertia matrix defined by 



Ji i 



This result is demonstrated in Woodhouse (16| 




y-j / A A A A \ 

R(o, 



Figure 2: Schematic diagram depicting the position vectors x and r in terms of a more 
appropriate reference, R(o, e l5 e 2 , e 3 ), which is embedded in the rigid body in such a way 
that the origin, o, and centre of mass coincide as do the basis vectors and the principal 
moments of inertia. 



Exploring this idea of a more appropriate reference further and restricting it to one, final 
choice allows still greater simplification; by specifying the basis vectors to coincide with 
the eigenvectors of the inertia matrix J(R) (i.e. the basis vectors are principal axes), 
equation (0) reduces to 



J x \UJi 
^33^3 



(3) 



Conservation of Angular Momentum in Terms of this More Appropriate Ref- 
erence 



Rewriting equation (|1|) in terms of this new reference leads firstly to 

+ wAq = E(r A /) (by the Coriolis theorem), 
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which in terms of the relation ([!]), is 

Ju ~dT + ^ 33 ~~ J22 ^ UJ2UJ ^ = S ( r A /) ' e i 
J ^~^ + (Ju ~ ^33)^3^1 = S(r A /) • e 2 

where the ( no contraction) are the principal moments of inertia defined by 

pr*dTL rb . 



f 

Ja, 



This description is in terms of a reference embedded within the rigid body in such a way 
that the origin and centre of mass coincide, as do the axes and principal moments of 
inertia. Failure to recognise this can lead to an incorrect formulation of the forces acting 
on the rigid body. 



3 A Free Surface Description Ignoring Surface Ten- 
sion 

It is assumed that the vertical position of the free surface can be written as a function 
of horizontal coordinates and time, that is, £3 = h(x 1 ,X2,t). The free surface can be 
described equivalently by the condition 

r)(x, t) = h(xi, x 2 , t) — £3 = 0. 

Requiring that the material derivative of this free surface to be zero amounts to writing 

-^ + Vh-[v 1 ,v 2 ]=v 3 . (4) 

A physical interpretation of this mathematics is that the free surface is an infinitesimally 
thin, perfectly elastic skin by which the fluid is contained. The free surface is in actuality 
not a distinct material entity. The approach is consequently slightly limited. For example, 
no convection of the free surface is allowed. Restricting the movement of surface nodes 
to vertical translations alone circumvents having to rewrite the above equations in terms 
of a completely general referential description. 

The variational equations are obtained in the usual manner. Multiplying by an arbitrary 
function w and integrating, the equations 

J w^dT + J wVh ■ [vi,V2]dT = J wv^dT 

are obtained. 

Equation (f|) is hyperbolic and one would therefore anticipate discontinuities or shocks 
in the solution. A standard Galerkin type formulation is not sufficient (see Hughes 0, 
Hughes, Mallet and Mizukami and Johnson ||) under such circumstances. 
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The Streamline-Upwind-Petrov-Galerkin (S.U.P.G.) Method with discontinuity captur- 
ing is one popular method resorted to when solving hyperbolic equations. This amounts 
to using 

W = W* + Ti(Vu>* • V) + T 2 (VW* ■ Vtangent), V tangent = Vh, 

where w* is the arbitrary function in the formulation, t\ and r 2 are constants described 
in Hughes, Mallet and Mizukami 0. 

Hermite elements are one more way in which to obtain smoother surfaces. Hermite 
elements are used for C 1 approximations which require non-zero first order derivatives 
at the nodes. 

Remark: Although surface tension is not considered significant in the free surfaces to 
be described in this work, the free surface equations to be used under such circumstances 
would be 

(<T1 - CT2) 'P^ 2GC ^ + VC 

where the left hand side is the total traction (the difference between that acting from 
within and that acting from without) acting on the interface, the last term is a tangential 
force contribution in instances where surface tension is non-constant, G contains the mean 
curvature 

d / (V77)i + (V77) 2 1 ] d / (V7y)i + (V7y) 2 ' 



dxi \ || V?7 || J dx 2 \ || V77 



and C is the surface tension. See Landau and Lifshitz Bird, Armstrong and 
Hassager |1 and Joseph and Renardy M in this regard. 



4 A Completely General Referential Description of 
an Incompressible, Newtonian Fluid 

The implications of free surfaces and fluid-rigid body interactions are a deforming fluid 
domain and the implementation of the majority of numerical time integration schemes 
then becomes problematic where a purely Eulerian description of the fluid motion has 
been used. The reason is that most numerical time integration schemes require succes- 
sive function evaluation at fixed spatial locations. Meshes rapidly snarl when purely 
Lagrangian descriptions are used. One alternative is to use a completely general referen- 
tial description. 

The transformation to the completely general reference only involves coordinates where 
used as spatial variables and the resultant description is therefore inertial in the same 
way as Lagrangian descriptions are. The conservation principles of mass and linear 
momentum remain the basis of the fluid model despite having resorted to a completely 
general reference (as opposed to a conventional Eulerian reference) for the description of 
the fluid. For the purposes of reading the equations which folow some slightly specialised 
notation will be introduced. 
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Notation 



In the equations that follow the ~ s merely indicate that quantities are functions of a 
deforming reference, x. The operators V and div are the referential counterparts of V 
and div, that is to say that 

~ d ~ d d d 

V = and dlv = ~^=- + + 

OX OXi OX2 OX3 

The symbol J is used to denote the determinant of the deformation gradient 

OX 

and the notation A : B is used to denote the matrix inner product A^B^. 

In terms of this newly established notation the conservation principles of linear momen- 
tum and mass may be written in primitive form as 



P 



+ VvF~\v-v ref ) 

ot 



J = pbJ + div P 



and 



Vv : F 1 



where P is the Piola-Kirchoff stress tensor of the first kind, P = &F * J. In terms of 
the constitutive relation for a Newtonian fluid, a = —pi + 2fxD. Thus 



-pi + p 



VvF 1 



VvF 1 



F-\j 



since 



2 



Vv+ Vv 



A comprehensive derivation of these equations is given in Childs |J. The derivation of 
a variational formulation is along similar lines to that for the conventional Navier-Stokes 
equations (the purely Eulerian description). 



For a fluid of constant density, the variational formulation 

~ -1 



/> 1 tv ■ ^-JdVt 



n 



dt 



P 



w ■ Vv 



F -(v- v re f) 



Jdtt 



P 



w ■ bJdVL 



pVw : F l Jdn 



+p I wPNdT 



2p I D{w) : T)(v)Jdn 



I qVv : F l dtt = 

is obtained, where q and w are respectively the arbitrary pressure and velocity of the 
variational formulation. Under the circumstances of a Dirichlet boundary condition the 
term involving the boundary integral may naturally be omitted. 
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5 Statement of the Total Problem in Dimensionless 
Form 

Suppose that length is to be measured in units of X and velocity in units of V (which is 
really parameterising time by T = y). Then 



x = xX, v = vV and t = t—. 



The equations governing the translation of the rigid body are accordingly divided through 
by j^r, the equations governing the rotation by p s X 3 V 2 , and the fluid equations by 

Using the symbol Ju( no SU m) to denote the ith, dimensionless principal moment of inertia 
of the rigid body, 



j _ Jii(no sum) 

Jii(no sum) 77k ; 

PsX b 



H to denote the transition matrix for a transition to a reference whose axes coincide with 
these principal moments of inertia, c to denote the centre of mass of the rigid body, T r b 
to denote the dimensionless surface of the rigid body, i to denote a dimensionless time, p/ 
and p s to denote the density of fluid and solid respectively, h to denote the dimensionless 
elevation of a free surface (h being a function of a reference prohibited from deforming in 
a lateral direction i.e. nodes belonging to the free surface move in the vertical only), v to 
denote the dimensionless velocity of the underlying fluid, F to denote the deformation 
gradient, J its determinant, v mesh to denote the dimensionless velocity of a mesh which 
is otherwise allowed to deform freely and Re to denote the Reynolds number, 



the combined, dimensionless, free surface-fluid-rigid body problem can then be stated as 
follows: find x rb , 0, v rb , cj, h and v (the dimensionless respective position, orientation, 
velocity and angular velocity of the rigid body, the elevation of the free surface and the 
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velocity field of the fluid), which satisfy 



Jll~j=- + W33 — ^22)^2^3 



EL 

Ps 



H f_ (x - c) A j-pJ + — Z)| n dT rb 



■ ei 



7 ^2 ^7 7 \ - - 

<J22—j^ + (yil — ^33 J^S 1 -^! 



£L 

Ps 



H j_ (x-c) A{-pI + —Djn dT rb 



e 2 



PL 

Ps 



H I (x - c) A \ -pi + — D \ n dT rb 



■ e 3 



Included 
> Rigid 
Body 



dv 



rb Pf 



dt mp s 



L{- pI+ l b }" df *» + v-> b 



dO 

~dt 

dt 



= OJ 



= v rb 



dh - Tt __ 1 

-r= + Vh- [v 1 ,v 2 \ = v 3 



Free 
Surface 



dv 



+ X7vF \v-v mesh ) 



X 

J = —bJ + divP 



-t 



Vv: F = 



where 



P = aF *J = 



(-pl+ — VvF 1 + (VvF y ) F l J 



Fluid 



subject to the "no slip" requirements 



v \v rb = v rb + U3 A (x - c) 



at fluid-rigid body interfaces and 
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at fixed, solid impermeable boundaries. Additional boundary conditions depend on the 
problem in question. 



6 Approximate Fluid Equation and its Implementa- 
tion 

Two topics which are relevant to the approximate fluid equation are dealing with the 
nonlinearity of the convective term and enforcing incompressibility in such a way as to 
eliminate pressure as a variable. Only against such a background can the approximate 
fluid equation be written. 



6.1 The Incompressibility Condition 

Two methods to enforce incompressibility and eliminate pressure as a variable are consid- 
ered relevant for review. The use of either a penalty method or an iterative, augmented 
Lagrangian method both have distinct, yet different advantages. 

The L.B.B. or B.B. condition (so named after Ladyzhenskaya, Babuska and Brezzi) must 
be taken into account in the implementation of either method. So-called "locking" or 
"chequerboard" modes can arise in the event of the L.B.B. condition not being satisfied. 



Work by authors such as Oden, Kikuchi and Song [0] can be referred to for a more 
in-depth treatment of the L.B.B. condition. 

A penalty method was used to eliminate pressure as a variable and enforce incompress- 
ibility in later examples, despite the acclaimed superiority of the iterative augmented 
Lagrangian approach (discussed shortly). The results obtained using the more efficient 
and more easily implemented penalty method compared favourably with those obtained 
by Simo and Armero ||15|| , who used the iterative augmented Lagrangian approach. 
In other examples the penalty method was used for reasons of standardisation (so that 
results could be compared against those cited in the literature). 



The Simple Penalty Method 



Penalty methods are closely related to the artificial compressibility method and provide 
a remarkably simple strategy for solving the equations. A fundamental assumption made 
in using the penalty method is that the incompressibility condition, 

div v = 0, 

can be replaced by one of small incompressibility, 

div v = —ep. 

where e is small. (It is tempting to substitute this expression for p in the momentum 
equations immediately, thereby eliminating the pressure from the equations altogether 
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at an analytic stage and prior to construction. This would circumvent the complicated 
construction which otherwise results. If such an approach were to be taken, however, the 
L.B.B. condition would not be satisfied and an appropriate underintegration rule would 
then be necessary to avoid an approximation which will lead to "locking" or "chequer 
board" modes.) 

The corresponding variational, small incompressibility is 

/ qdxvvdCl = — e / qpdVt 
Jn Jn 

where q is the arbitrary function of the variational formulation. A well known shortcoming 
of penalty methods is the severe ill-conditioning which occurs as the penalty parameter 
is decreased, making the simulation of exact incompressibility practically impossible. As 
the theoretical limit e = 0, an infinite condition number for the discrete algebraic sys- 
tem is obtained. Despite this shortcoming a straight forward penalty method provides 
a remarkably easy-to-implement means of enforcing incompressibility and eliminating 
pressure as a variable, furthermore, there was no discernable difference between results 
obtained using the penalty method and those obtained using the iterative augmented La- 
grangian approach for the driven cavity flow problem (compare the results in Subsection 
Owith those obtained by Simo and Armero ]T5|). 



An Iterative Augmented Lagrangian Approach for Enforcing Exact Incom- 
pressibility 

There is considerable merit in the iterative augmented Lagrangian approach from the 
point of view of realising theory. Incompressibility can be exactly enforced (to within a 
specified tolerance) for finite values of the penalty parameter using an iterative augmented 
Lagrangian approach. This is in contrast to conventional penalty formulations where 
there is a trade-off between enforcing exact incompressibility and obtaining a severely 
ill-conditioned algebraic system. 

A fundamental assumption made in using the iterative augmented Lagrangian method 
is that the incompressibility condition can be replaced by one of iterative, limiting-case 
incompressibility, 



div?; = — 




for e small. 

The method of augmented Lagrangians provides an effective means to circumvent the 
shortcomings of the simple penalty method and enables the divergence free constraint to 
be exactly enforced for finite values of the penalty parameter. 



6.2 A New, Linear Scheme for the Approximation of the Con- 
vective Term 

The following theorem is the first of two proposed numerical improvements (originally 
presented in Childs and Reddy Linearising with a guess obtained by extrapolating 
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through solutions from the previous two time steps leads to second order accuracy. This 
is an improvement on the conventional method by an order of magnitude. 



Theorem 1 The linearised terms, (2v \t —v \t-At) ■ Vv \t+At and v \t+At -V(2v \t 
—v \t-At), are second order accurate (have error 0(At 2 ) ) approximations of the nonlinear 
term (v ■ Vu) \ t +At- 



Proof: 



V | t+At 



, dv 
v t +At — 
ot 



0{At 2 



(by Taylor series) 



v |t +At 



V |t -V \t-At 

At 



+ O(At) + 0(At 2 ) (using a backward 



2v \ t -v |t-At +0(At 2 



difference) 



{v ■ Vv) |t+At = [2v \ t -v |t-Ai +0(At 2 

= [2v \ t -v |t-At] • | t+At +0(At 



■ VV |t+At 
2\ 



The above linearisation schemes are an improvement on the conventional v \ t • Vf |t+At H 
or v |t+At "Vf \t linearisation schemes by an order of magnitude. 

Further refinement of this second order accurate solution can usually be accomplished 
by the process of Picard iteration, if necessary - providing the initial guess lies within a 
radius of convergence. Picard iteration amounts to "linearising" the nonlinear term with 
the solution obtained during the previous iterate. 



6.3 Approximate Fluid Equation 

An updated approach and the choice of a referential configuration which coincides with 
the spatial configuration at the instant within each time step about which the equations 
are to be evaluated, facilitates the elimination of the deformation gradient from the re- 
sulting fluid scheme completely (under such conditions the deformation gradient becomes 
identity - see Childs ||). 

A backward difference is used to approximate the time derivative in accordance with the 
findings of Childs 0. The finite element method is used for the spatial (referential 
"space" - the reference deforms) discretisation and a Q2-P1 element pair is the chosen 
basis. A penalty method is used to eliminate pressure as a variable. The approximate 

^Favoured in terms of both rate and radius of convergence by Cuvelier, Segal and van Steen- 

HOVEN §. 
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equations, to be solved for the nodal velocities, (pertaining to element e), are conse- 
quently 



A { -h 



BkiB k jJ e dVl + / B ki G k jiB lm J e dVt{V r 



e linearisation pre mes/i\ 



+ 



i -1 



Lp m LpiJ e dVL 



. l PmG k j k J e dVL 1^ (fiG n i n J e dVt 



= A B ki B k3 J e dCl + ± ^S fei S fcj J e df2(^ e | t _ At )}, 

where A is the element assembly operator, i? is the total number of elements, e, into 
which the domain has been subdivided, Cl is the master element domain, At is the length 
of the time step, the are the basis functions, {£} is the master element coordinate 

system, 



B(0 



0! ... n 

0i ... n 
0! ... 0, 



¥»(€) = [l,xi(€) ) a: a (£) ) X3(£)]I > 



f dec 1 

(£) , J e = det | — > for element e, 







n is the number of nodes on each element, /i is viscosity, p is density, X, V are a length 
scale and a velocity scale respectively, V e hneansaUon j s the new second order accurate 
linearisation (given in Subsection |6.2|) , 



'e linearisation 



2V e \ t -V e 



It- At, 



V e mesh is the mesh velocity, e is the penalty parameter (e = 10 6 ) and g is the gravita- 
tional acceleration. 

To recover the pressure on each element once the velocity solution has been obtained, 



V 



(p m ipiJ e dtt 



fmG kjk J e dfl V e pi. 



(5) 



6.4 A Guideline for Artificial Reynolds Number Adjustment 

A straight forward Picard iteration process (in which the convective term is linearised) 
may not be sufficient at higher Reynolds numbers where the nonlinear term becomes 

tThe Q 2 -Pi element pair was shown to satisfy the L.B.B. condition in the context of rectangular 
elements. It is important to note in this regard that a linear function mapped from the master element 
using a Q2 mapping will no longer be Pi for non-rectangular elements. 
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more significant and the radius of convergence is smaller. What one has to do under 
such circumstances is implement a continuation technique. That is, lower the Reynolds 
number to a level where convergence is attainable, then raise it incrementally, obtaining 
progressively better solutions at each stage. Choosing these various artificial Reynolds 
number levels can, however, turn out to be somewhat problematic. What follows was 
devised as a guideline in making this choice. Suppose the nonlinear and linear parts 
of the equations are expressed in terms of two matrices, A and B, respectively. The 
equations to be solved are then 

A n kjV k Vj + B nj Vj + c n = 0. 



Theorem 2 A sufficient condition for the iterative linearisation method to converge is 



D \\ p < 1 



where 

Dji = (Ankjvl. — B n j) A n [ m v l m , 



| . ||p denotes the generalised p norm and the superscript i denotes the number of the 
iteration from which a given solution was obtained ('Childs and Reddy j^j). 

PROOF: Let e J_1 denote the difference between the solutions obtained from the ith and 
(i — l)th iterations. That is, 

e*- 1 = v i - 

The system of equations solved at the (i — l)th Picard iterate was accordingly 

AnkjV^vir 1 + B nj v}~ 1 + c n = 
A nkj (vl- 1 -ei- 2 )vy 1 +B nj v^ l + c n = 0. (6) 

The system of equations to be solved at the ith Picard iterate is 

AnkjV^v) + B nj v) + c n = 
A nkj vl 1 (v}- l + ey 1 ) + B nj (v}- l + ey 1 ) + c n = 0. (7) 

Equating (§) with (0), 

A — 1 I p A— 1 A A — 2A—1 

s^nkjUk t j "r n nj t j — ~ n -nkj t k U j 

e)- 1 = -(A^vl' + B^A^el^- 1 

J- 



D - 1 

S— 1 II ^ 1 1 n 1 1 II -i— 2 



\\p — ii i ip 1 1 *" i ip ■ 



A sufficient requirement for convergence is therefore 



€ I Ip 11^ 1 1 p 1 1 1 1 p < ^ 
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Corollary 1 A sufficient condition for the iterative linearisation method to converge 
is 

max | {A^v]^ 1 - B nj )~ l A nlm v^ 1 |< 1 
i J 

(the maximum absolute column sum < 1), where the superscript i denotes the number of 
the iteration from which that particular solution was obtained. 

PROOF: In the particular instance of the 1 vector norm, 

|| D ||i= max £™ | (A nk jvj~ 1 - B nj )~ l A^v^ 1 \ 
i J 

(The natural matrix norm induced by the 1 vector norm is the maximum absolute column 
sum.) 

Remark: This guideline serves only to examine local convergence, and one still has to 
hunt for a suitable Reynolds number by trial and error. It does, however, is provide an 
exact quantitative measure of this local convergence and, depending on the solver used, 
can improve efficiency 



7 A Predictor— Corrector or "Linearised" Approach 

The combined free surface-fluid-rigid body problem is highly nonlinear. An important 
decision had to be taken on whether to solve the respective sub-problems simultaneously, 
using a Newton-Raphson method, or separately using a predictor-corrector (iterative) 
approach. The latter approach may be justified in terms of a Picard-type linearisation 
process (similar to that discussed in Section |6^ when dealing with the convective term). 

The advantages of a predictor-corrector approach over solving simultaneously (the ap- 
proach taken by many) are expected to be a larger radius of convergence and hence a 
more robust scheme, substantially less code and easy, separate testing of sub-problems. 

The disadvantages of a predictor-corrector approach over solving simultaneously are ex- 
pected to be linear, as opposed to quadratic, rates of convergence. The expected slow 
rate of convergence could possibly be overcome to a degree by taking good initial guesses 
based on previous solutions, using linear extrapolation for example. 

7.1 The Resulting Algorithm 

Figure |3] outlines the algorithm which results from the use of such a linearised approach. 
The solution values of the fluid sub-problem are used in solving the free surface sub- 
problem, working out the traction on the rigid body and consequently in solving the rigid 
body sub-problem. The mesh is then adjusted accordingly and, if not the first iteration, 
the convergence is examined. If only the first iteration, alternatively, if the degree of 
convergence is found to be unsatisfactory, the process is repeated. Once this iteration 
process has converged, a new time step is commenced. 
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The 0(At 4 ) Runge-Kutta-Fehlberg and Adams-Moulton approximations used in the 
rigid body sub-problem are possibly excessive in comparison to the O(At) finite difference 
schemes implemented in the fluid and free surface problems. 



Flow of Variables 
and Sequence 

Flow of Variables 
Only 

Sequence Only 



Set "Start Up" 
Values/Guesses, 



Basic Mesh 



Tailor Mesh 







Workout Skyline 
Storage 







Set "Previous" 
Values = "Present" 
Values 



Set Dynamic Bou- 
ndary Conditions, 
t = t + dt etc. 



Solve Fluid 
Equations 



• ^ ^* Traction and Trac- 
tional Torque on 
Rigid Body 



Solve Free Surface 
Equations 



Solve Equations of 
Motion 



"Previous iteration 
values = "present 
iteration" values 



Write Results Out , 




Figure 3: Combined Structure-Flow Chart Outlining the Predictor-Corrector Algorithm 
for the Fluid-Rigid-Body-Free Surface Problem. 
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8 Automatic Mesh Generation About Included Bod- 
ies of Any Shape 

The finite element mesh was automatically generated and adjusted about the included 
rigid body in what is possibly a slightly novel fashion. A small region of mesh immediately 
adjacent to the included rigid body was repeatedly remapped to cope with the changing 
orientation, the remainder was squashed/stretched according to the translation. 



Local Distortion About an Included Rigid Body 

The local distortion about included rigid bodies is obtained by mapping four square 
chunks of rectangular mesh to the four wedge-shaped domains depicted in Figure f| using 
finite element mappings. A square region of mesh centered on, and including the rigid 
body, is first removed. Each of the depicted wedge-shaped regions is then demarcated 
by as many points as there are nodes in an element i.e. each wedge shaped-region is set 
up as a massive element. 

Demarcation: Demarcation is accomplished by first locating the intersection of the lines 
which bisect corners and edges of the square frame, with the surface of the rigid body 
using Newton's method. The remainder of this demarcation process needs no explanation. 




(-1. i) (h i) 




(-1, -1) (I, -1) 

Figure 4: The Local Distortion is Obtained by Mapping Square Chunks of Rectangular 
Mesh Using Finite Element Mappings. 
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Mapping: Chunks of uniform mesh, which have identical extremities to those of the 
master element, are then mapped into the newly-demarcated, wedge-shaped regions 
using finite element mappings. The devised method maps chunks of uniform mesh into 
wedge-shaped regions surrounding the included rigid body in exactly the same manner 
as points in the master element domain are, in theory, mapped into individual mesh 
elements. In the finite element mapping, 

nNode 

i=i 

symbols which previously denoted the master element, the master element variables, the 
corresponding position in the domain, the position of node i and element e now denote 
the following: Cl - the chunk of uniform mesh, £ - the position of a constituent node in 
the chunk, x - the corresponding position in the wedge, Xi - the position of the ith point 
demarcating the wedge and Q e - the wedge. 

Interface Adjustment: Further fine adjustment of nodes intended to delineate the 
surface of the rigid body is accomplished by moving them along a line between node 
and centre, to the rigid body surface using Newton's method, thereby accomodating 
additional shape complexity. 




Figure 5: Automatic Mesh Generation about a Rigid Body which is Simultaneously 
Rotating and Translating. 

The method described was found to be remarkably practical, simple and effective with 
maximum angles never exceeding j radians at element vertices within the mesh. 
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Translational Mesh Deformation 

The mesh outside the "box" (the box containing the 4 "wedges" enclosing the rigid 
body) is squashed/stretched according to the requirements of the translation (the nodes 
are translated in the same direction by a factor inversely proportional to their distance 
from the box). 




Figure 6: Automatic Mesh Generation and Adjustment about a Diamond which is Si- 
multaneously Rotating and Translating. 



9 Some Test Examples 

The final model was tested in the context of a driven cavity flow, a driven cavity flow 
with various, included rigid bodies, a die-swell problem and a Stokes, second-order wave. 
The idea was to use as simple as possible flows at extremely low Reynolds numbers, so as 
to generate the smoothest possible flows in which any resulting rigid body motion could 
be expected to be highly predictable, alternatively examples for which there exist well 
established results. 



9.1 Example 1: Driven Cavity Flow 



A variety of bench-mark tests for viscous, incompressible flows are cited in the literature. 
The driven cavity flows of Carey and Krishnan [g] and Simo and Armero |L5| are two 
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such examples. The so-called "leaky lid" driven cavity flow of Simo et. al. was found to 
be the more sensitive of the two when testing the performance of schemes approximating 
the Navier-Stokes equations. In addition, Simo et. al. tested the scheme for the 
Navier-Stokes equations rigorously using the leaky lid problem. 

These results could be used as a comparison in testing the performance of the completely 
general referential description derived in Childs || as well as to compare the perfor- 
mance of the standard, penalty method (used here) with that of the iterative, augmented 
Lagrangian approach (used by Simo et. al.). The main objective of these tests was 
nonetheless to generate and study as smooth as possible a flow into which a somewhat 
inconsequential (flow-wise inconsequential) rigid body could later be introduced (see the 
forthcoming example). 

Remark: Both the above mentioned flows are steady in the context of the traditional 
Eulerian description. They are, however, "unsteady" in the context of a deforming mesh, 
a factor enhancing their potential as tests. 

The problem is that of a square, two-dimensional "pot" whose lid is moved across the top 
at a rate equal to its diameter for a Reynolds number of unity. The boundary conditions 
are accordingly "no slip" on container walls and a horizontal flow of unity across the top 
(as depicted in Figure 0). The problem also bears a certain resemblance to an idealised 
pothole in a river bed. A mesh consisting of one hundred and forty four elements was 
deformed at differing rates (also depicted in Figure |7|). 




6 The Leaky Lid Problem 1 The uniform mesh ( 144 elements). Navier-Stokes pressure (control). 



Figure 7: The Problem, the Mesh and the Pressures Obtained Using the Conventional 
Eulerian Equations. 

Results 



The results in Figures ^ and |10| obtained while deforming the 144 element mesh de- 
picted, are in agreement with those obtained using the analogous Navier-Stokes algorithm 
on a fixed uniform mesh (the control results in Figure |7| and the relevant velocity profiles 
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in Figure [T^ on page ^3|). They are also in agreement with the results of Simo and 
Armero obtained using the Navier-Stokes equations (the conventional Eulerian 
description) on a fixed, uniform mesh consisting of 3200 triangular elements. 




The uniform mesh ( 144 elements) 
at the start of the test time step. 
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The distorted mesh ( 144 elements) 
at the end of the test time step. 



Pressure for a 5% mesh 
ssion. 



compre- 



Figure 8: Pressures Obtained Using the Completely General Referential Equation. 



The fixed, uniform mesh, with which the control results were obtained, consisted of 144 
9-noded, quadrilateral elements. In the first test all, but the 25 upper boundary nodes, 
were compressed into the lower 95 % of the problem domain after which the mesh was 
restored to a uniform state. 




The distorted mesh ( 144 elements) The uniform mesh ( 144 elements) Pressure for a 5% mesh decompre- 
at the start of the next time step. at the end of the next time step, ssion. 



Figure 9: Pressures Obtained Using the Completely General Referential Equation. 



The mesh was successively compressed and decompressed over two time steps, each of 
length 0.05. Mesh velocities were therefore of the order of 10 times greater than the flow 
velocities being modelled. 
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Figure 10: In this test part of the mesh was successively compressed and decompressed 
by 5 % over two time steps of length 0.05. 
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In the second test (Figs. [0], [12| and [13]) all but 25 boundary nodes were compressed into 
the lower three quarters of the domain after which the mesh was restored to a uniform 
state. 




The umfojin mesh £144elemenlEj The dilcj led mesh (144 element J FjesEiue fci a 25'SjjneEh can pe- 
31 ihe Etui cf ihe lesl lime dep. al ihe end of ihe 1ce1 time Elep. EEion. 



Figure 11: Pressures Obtained Using the Completely General Referential Equation. 



In this example the mesh was successively compressed and decompressed by 25 % over 
two time steps of length 0.05. Mesh velocities were therefore of the order of 50 times 
greater than the flow velocities being modelled. 





The diElQjIed mesh i 144 eLeme rrlE J The u nifcim me Eh ( 144 eLeme nlE 3 
al ihe Elajl <if ihe nez1 limeElep. al iheendcf ihe nez1 limeElep. 



1 2 £ E I £ 3 £ 3 n. 
PjeEEUje foj a 25* mesh decam- 

pjeEElCfl. 



Figure 12: Pressures Obtained Using the Completely General Referential Equation. 



The error is attributed to a combination of a bad mesh (mostly) and relatively high 
off-scale mesh velocities. It should be remembered that a bad mesh is a bad mesh 
nonetheless. 
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Figure 13: In this test part of the mesh was successively compressed and decompressed 
by 25 % over two time steps of length 0.05. 
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Remark: A further fact illustrated by the results of this test is that, of the pressure 
and velocity solutions, the pressure solution is the more sensitive of the two. This is 
to be expected since the pressure is recovered from what is fundamentally a difference 
in something like the tenth decimal place of the velocity solution (see equation (|[)) 
depending on the penalty parameter. 



9.2 Example 2: "Pebble in a Pothole" 

In this example rigid bodies of varying mass and moments of inertia were released from 
rest in a flow dictated by the same boundary conditions as the driven cavity flow of 
the previous example. One would expect a "die bead" (a small rigid body of neutral 
bouyancy) to move in tandem with the fluid soon after its release from rest. One might 
also expect a clockwise spinning to be induced by concentrating the mass closer to the 
centre i.e. lowering the moment of inertia. 

Various rigid bodies were introduced to the cavity driven flow problem described in 



Subsection |9J], in the absence of a body force. Although the algorithm allows for the 
inclusion of a rigid body of any shape, one which could be expected to behave fairly 
predictably was selected for the purposes of the test. The results presented shortly 
involve the elliptical rod-shaped body 

| + ft = 0-025 2 , 

whose major axis is 0.1. 



Results 



The results which follow were in agreement with the notion that the die bead would move 
in tandem with the fluid soon after its release from rest. In the succesive trajectories in 



Figure [14] the mass was successively concentrated closer to the centre (a lower moment 



of inertia was used). A clockwise spin was then induced. 

The last trajectory in Figure [14] was obtained using a Reynolds number of unity in order 
to generate a smooth, predictable flow as close as possible to the driven cavity flow 



of Subsection [14]. Although the scaling is not immediately reminiscent of any real life 
problem, the example serves to further verify the potential of the methods described 
qualitatively. 

The results for the included rigid body tests were extremely encouraging. This is es- 
pecially so when it is considered that the Reynolds numbers implied by small, included 
rigid bodies are low and the model being developed is eventually intended to elucidate 
problems of a sedimentological nature. 
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Figure 14: The trajectories of various included rigid bodies released from rest at the 



centre of the driven cavity flow described in Subsection 9.1. Top Left: Re = 0.025 



m = 251.3, J 33 = 314.2 and t = 3.6 sees. Top Right: Re = 0.025, m = 251.3, J 33 = 1.0 
and t = 4.0 sees. Bottom Left: Re = 0.025, m = 251.3, J 33 = 0.1 and t = 3.6 sees. 
Bottom Right: Re — 1, m — 1, moment of inertia (scaled) = 0.1 and t = 2.0 sees. 



9.3 Example 3: Die Swell Problems 

The axis-symmetric die swell (or fluid jet) problem is a free surface problem well doc- 
umented in the literature. The exact details of the problem description always differ to 
varying degrees, however, the central theme basic to all involves the extrusion of a fluid 
with initial parabolic flow profile from the end of a short nozzle. The various boundaries 



in terms of which this particular version is specified are depicted in Figure [15. 
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Figure 15: The Die Swell Problem 



At the free surface, T: 



At the nozzle wall, I\ 
At the inlet, Tf. 



At the outlet, r o : 



h 

dh 
dxi 



At the symmetry axis, T s : 



x 1= 



x\=L 



V 



V2 



an ■ n 



X 




Vl = 



2 V X 2 



vt = 



an ■ t 



v ■ n = 



3 fi ( x 2 2 " 



o 

(= PN-N) 

(= v-T) 

(= PN-T) 

(= v-N) 
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The Symmetry Axis, T s : Simultaneous implementation of both T s boundary condi- 
tions amounts to leaving the term 



w ■ PN dT, 



out of the fluid equations and specifying v ■ N, the component of velocity normal to 
the boundary. Because the boundary condition is Dirichlet in the normal direction, 
variational formulation is consequently not required in that direction. This allows w in 



/ w ■ PN dT s = f (PN ■ N)(w ■ N) dT s + f (PN ■ T)(w ■ T) dT s (8) 



to be chosen arbitrarily, but in such a way that w ■ N — 0. 

The Outlet, r o : Implementing both r o boundary conditions simultaneously amounts 
to omitting the term 

w ■ PN dT 

To 



and specifying v ■ T. This implementation can be justified in a likewise fashion to that 
given for the boundary condition at the symmetry axis. 



Results 



Examples of die swell problems in which X = 1, L Q = 3.5 and L = 20 for various 
kinematic viscosities and/or velocity scales were examined using a 150 (50 x 3) element 
mesh. The predicted die swell ratios agree well with those obtained by other, more 
specific methods. 



Compare [] the result of Omodei [TJ], Figure 5 (on page 87) with the corresponding 



Re = 1 result depicted in Figure 16, for example. Good agreement is obtained between 



Omodei's Figure 7 (on page 88) result and the corresponding Re = 10 result depicted 



in Figure [17], likewise between his Figure 8 (on page 89) result and the corresponding 



Re = 18 result depicted in Figure [TBI. These results are further substantiated by those of 



Kruyt, Cuvelier, Segal and Van Der Zanden |pT| , and those of Engelman and 
Dupret quoted in Kruyt et al.. 



1 Care needs to be taken in matching inlet velocity profiles as well as scales when comparing results. 
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Figure 16: Die swell ratios predicted for various Reynolds numbers using an inlet velocity 
profile of vi = 



■(1 — x|) and the methods described. 
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Figure 17: Die swell ratios predicted for various Reynolds numbers using an inlet velocity 



profile of vi 



Re 3 
10 2 



(1 — x%) and the methods described. 
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Figure 18: Die swell ratios predicted for various Reynolds numbers using an inlet velocity 
profile of V\ — 



2l| 2 (1 — x|) and the methods described. 



9.4 Example 4: A Stokes Second Order Wave 



In this test the velocity profile and surface elevation predicted by Stokes second order 
wave theory were used as boundary conditions for flow and free surface subproblems 
respectively. The development of Stokes second order theory is the same as that of 
Airy wave (Stokes first order) theory, the exception being that terms as far as second 
order (O 



H 2 



x 2 



in the power series expansion of the velocity potential are used. Wave 
troughs are consequently flatter and the crests are sharper than those predicted using 
the corresponding first order theory (see KOUTITAS [TO] in this regard). 



Results 



The method converged rapidly to the Stokes second order surface (Figure [19| on page 
32|) in the initial stages, despite the use of some fairly preposterous free surface starting 
guesses. This was found to be so for a range of wave lengths, amplitudes and domains 
experimented with. The problem was, however, not attempted with the same serious- 
ness as previous examples and the mesh was poor (there being only three elements in 
the vertical extent of the mesh). Problems with wave propagation were consequently 
experienced as time progressed. 
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10 Conclusions 



The completely general referential description proposed in Childs || is varified by the 
numerical experiments of Section |9|, as was the predictor-corrector approach proposed 
in Section [7[ Although the resources at the authors disposal limited the tests to two 
dimensional problems, this is not expected to be of any consequence. 

The linearised terms (2v \t —v \t-At) • Vi> \t+At and v \t+At - V(2i; | t — v \t-At) are second 
order accurate approximations of the convective term. The formula given in Theorem ^ 
provides an exact quantitative measure of local convergence for a given Reynolds number 
when applying a continuation technique and, depending on the solver used, can improve 
efficiency 

The more efficient and more easily implemented penalty method was found to produce 
comparable results to the iterative augmented Lagrangian approach (through which in- 
compressibility can be exactly enforced) in the simple test examples performed. With 
regard to the L.B.B. condition, it is important to note that a linear basis function mapped 
from the master element using a Q2 mapping will no longer be Pi for non-rectangular el- 
ements (the Q2-P1 element pair was shown to satisfy the L.B.B. condition in the context 
of rectangular elements). 

The holistic method to automatically generate meshes about rigid bodies included within 
the fluid, using finite element mappings, was found to be remarkably practical, simple 
and effective with maximum angles within the mesh never exceeding | radians (outlined 
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in Section 
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